Stochastic Hard-Sphere Dynamics for Hydrodynamics of Non-Ideal Fluids 
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A novel stochastic fluid model is proposed with non-ideal structure factor consistent with com- 
pressibility, and adjustable transport coefficients. This Stochastic Hard Sphere Dynamics (SHSD) 
algorithm is a modification of the Direct Simulation Monte Carlo (DSMC) algorithm and has several 
computational advantages over event-driven hard-sphere molecular dynamics. Surprisingly, SHSD 
results in an equation of state and pair correlation function identical to that of a deterministic 
Hamiltonian system of penetrable spheres interacting with linear core pair potentials. The fluctuat- 
ing hydrodynamic behavior of the SHSD fluid is verified for the Brownian motion of a nano-particle 
suspended in a compressible solvent. 
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With the increased interest in nano- and micro-fluidics, 
it has become necessary to develop tools for hydrodynamic 
calculations at the atomistic scale 111 Of particular inter- 
est is the modeling of flexible polymers in a flowing solvent 
for both biological (e.g., cell membranes) and engineering 
(e.g., micro-channel DNA arrays) applications. Typically 
the polymer chains are modeled using Molecular Dynamics 
(MD). For many applications, a realistic representation of 
the solvent and bidirectional coupling between the flow and 
the polymer motion is needed, for example, in the model- 
ing of turbulent drag reduction. Previously, we introduced 
the Stochastic Event-Driven MD (SEDMD) algorithm that 
uses Direct Simulation Monte Carlo (DSMC) for the solvent 
coupled to deterministic EDMD for the polymer chain 
However, DSMC is limited to perfect gases. Efforts have 
been undertaken to develop solvents that have a non-ideal 
EOS, and that also have greater computational efficiency 
than brute-force molecular dynamics. Examples include 
the Lattice-BoltzmannjXB) method Dissipative Parti- 
cle Dynamics (DPD) [5(, and Multi-Particle Collision Dy- 
namics (MPCD) Qjeach of which has its own advantages 
and disadvantages [1]. The Stochastic Hard Sphere Dy- 
namics (SHSD) algorithm described in this Letter is based 
on successive stochastic collisions of variable hard-sphere 
diameters and is thcrmodynamically consistent (i.e., the 
direct calculation of compressibility from density fluctua- 
tions agrees with the density derivative of pressure) . SHSD 
modifies previous algorithms for solving the Enskog kinetic 
equation 0, Q while maintaining good efficiency. 

In the SHSD algorithm randomly chosen pairs of ap- 
proaching particles that lie less than a given diameter of 
each other undergo collisions as if they were hard spheres 
of diameter equal to their actual separation. The SHSD 
fluid is shown to be non-ideal, with structure and equation 
of state equivalent to that of a fluid mixture where spheres 
effectively interact with a repulsive linear core pairwise po- 
tential. We theoretically demonstrate this correspondence 
at low densities. Remarkably, we numerically find that this 
effective interaction potential, similar to the quadratic core 
potential used in many DPD variants, is valid at all densi- 
ties. Therefore, the SHSD fluid, as DPD, is intrinsically 
thermodynamically-consistent, while non-ideal MPCD is 
only numerically thermodynamically-consistent for tuned 
choices of the parameters jy, 0] . 



As an algorithm, SHSD is similar in nature to DPD 
and has a similar computational complexity. In DPD, mo- 
mentum is also stochastically exchanged between particles 
closer than a given distance. The essential difference is 
that DPD has a continuous-time formulation (a system of 
stochastic ODEs), where as the SHSD dynamics is discon- 
tinuous in time. This is similar to the difference between 
MD for continuous potentials and discontinuous potentials. 
Just as DSMC is a stochastic alternative to hard-sphere MD 
for low-density gases, SHSD is a stochastic modification of 
hard-sphere MD for dense gases. On the other hand, DPD 
is a modification of MD for smooth potentials to allow for 
larger time-steps and a hydrodynamically-consistent ther- 
mostat. 

The SHSD algorithm is not as efficient as DSMC at a 
comparable collision rate. However, when low compress- 
ibility is desired, SHSD is several times faster than EDMD 
for hard spheres, the fastest available deterministic alter- 
native. Low compressibility, for example, is desirable so 
that flows are kept subsonic even for high Reynolds number 
flows. Furthermore, SHSD has several important advan- 
tages over EDMD, in addition to its simplicity: (1) SHSD 
has several controllable parameters that can be used to 
change the transport coefficients and compressibility, while 
EDMD only has density; (2) SHSD is time-driven rather 
than event-driven thus allowing for easy parallelization; (3) 
SHSD can be more easily coupled to continuum hydrody- 
namic solvers, just like ideal-gas DSMC [l(|. Strongly- 
structured particle systems, such as fluids with strong in- 
terparticle repulsion (e.g., hard spheres), are more difficult 



to couple to hydrodynamic solvers [ll| than ideal fluids, 
such as MPCD or DSMC, or weakly-structured fluids, such 
as DPD or SHSD fluids. 



The standard DSMC [12| algorithm starts with a time 
step where particles are propagated advectively, r i = + 
ViAt, and sorted into a grid of cells. Then, a certain 
number N co u ~ T SC N C (N C — l)At of stochastic collisions 
are executed between pairs of particles randomly chosen 
from the N c particles inside the cell. The conservative 
stochastic collisions exchange momentum and energy be- 
tween two particles i and j that is not correlated with the 
actual positions of the particles. Typically the probability 
of collision is made proportional to the magnitude of the 
relative velocity v r — |vy| by using a conventional rejec- 
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tion procedure. DSMC, unlike MD, is not microscopically 
isotropic and does not conserve angular momentum, lead- 
ing to an anisotropic stress tensor. To avoid such grid ar- 
tifacts, all collision partners within a collision diameter D 
must be considered even if they are in neighboring cells, 
and, if angular momentum conservation is required, only 
radial momentum should be exchanged in collisions as for 
hard spheres. This grid-free variant will be called Isotropic 
DSMC (I-DSMC). The cost is that is the computational 
efficiency is reduced by a factor of 2 — 3 due to the need 
to perform neighbor searches. Note that a pairwise Ander- 
son thermostat proposed within the context of MD/DPD 
in Ref. [3] essentially adds (thermostated) I-DSMC colli- 
sions to ordinary MD and has very similar computational 
behavior. As in I-DSMC, in SHSD we consider particles 
in neighboring cells as collision partners in order to en- 
sure isotropy of the collisional (non-ideal) component of 
the pressure tensor. 

The virial (Avy • Ar.y) vanishes in I-DSMC giving an 
ideal-gas pressure. In order to introduce a non-trivial equa- 
tion of state it is necessary to either give an additional 
displacement to the particles that is parallel to Avy, or 
to bias the momentum exchange Av.y to be (statistically) 
aligned to Ar ^ . The former approach has already been in- 
vestigated in the Consistent Boltzmann Algorithm (CBA) 
however, CBA is not thermodynamically consistent 
since it modifies the compressibility without affecting the 
density fluctuations (i.e., the structure of the fluid is still 
that of a perfect gas). A fully consistent approach is to 
require that the particles collide as if they are elastic hard 
spheres of diameter equal to the distance between them at 
the time of the collision. Such collisions produce a posi- 
tive virial only if the particles are approaching each other, 
v n — — Vy- • f ij > 0, therefore, we reject collisions among 
particles that are moving apart. Furthermore, as for hard 
spheres, it is necessary to collide pairs with probability that 
is linear in v n , which requires a further increase of the re- 
jection rate and thus decrease of the efficiency. Without 
rejection based on v n or v r , fluctuations of the local tem- 
perature T c would not be consistently coupled to the local 
pressure p c ~ (Av^ • Ar^) ~ T sc \/T~ c because p c would be 
~ v^c instead of the necessary p c ~ T c . For DSMC the col- 
lisional rules can be manipulated arbitrarily to obtain the 
desired transport coefficients, however, for non-ideal fluids 
thermodynamic requirements eliminate some of the free- 
dom. This important observation has not been taken into 
account in other algorithms that randomize hard-sphere 
MD Note that one can in fact add I-DSMC collisions 
to SHSD in order to tune the viscosity without affecting 
the compressibility. 

For sufficiently small time steps, the SHSD fluid can be 
considered as a simple modification of the standard hard- 
sphere fluid. Particles move ballistically in-between colli- 
sions. When two particles i and j are less than a diameter 
apart, ry < D, there is a probability rate (3%/ D)v n @(v n ) 
for them to collide as if they were elastic hard spheres with 
a variable diameter D$ — rij. Here 9 is the Heaviside 
function, and x is a dimensionless parameter determining 



the collision frequency. The prefactor 3/D has been chosen 
so that for an ideal gas the average collisional rate would 
be x times larger than that of a low-density hard-sphere 
gas with density (volume fraction) <fi = ttND 3 /(6V). 

In order to understand properties of the SHSD fluid as 
a function of <f> and Xj we consider the equilibrium pair 
correlation function gi at low densities, where correlations 
higher than pairwise can be ignored. We consider the cloud 
of point walkers ij representing the N(N — l)/2 pairs of 
particles, each at position r = — r, and with veloc- 
ity v = Vj — Vj. At equilibrium, the distribution of the 
point walkers in phase space will be /(v, r) = f(v r ,r) ~ 
g%{r) exp(— mv^/AkT). Inside the core r < D this distri- 
bution of pair walkers satisfies a kinetic equation 



df 
dt 



df 
or 



= v n Tof, 



where Tq = 3%/Z) is the collision frequency. At equilib- 
rium, df /dt = and v n cancels, consistent with choos- 
ing collision probability linear in \v n \. Thus dgi/dx = 
3x<?2@(l — x), with solution g2{x) — exp [3%(x — 1)] for 
x < 1 and g2(x) = 1 for x > 1, where x = r/D. In- 
deed, numerical experiments confirmed that at sufficiently 
low densities the equilibrium gi for the SHSD fluid has 
this exponential form inside the collision core. This low 
density result is equivalent to g% = exp[— U(r)/kT], where 
U(r)/kT = 3%(1 — x)Q(l — x) is an effective linear core 
pair potential similar to the quadratic core potential used 
in DPD. Remarkably, it was found numerically that this 
repulsive potential can predict exactly g% (x) at all liquid 
densities. Figure [T] shows a comparison between the pair 
correlation function of the SHSD fluid on one hand, and 
a Monte Carlo calculation using the linear core pair po- 
tential on the other, at several densities. Also shown is a 
numerical solution to the hyper-netted chain (HNC) inte- 
gral equations for the linear core system, inspired by its 



success for the Gaussian core model [16|. The excellent 
agreement at all densities permits the use of the HNC re- 
sult in practical applications, notably the calculation of the 
transport coefficients. 
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Figure 1: (Color online) Equilibrium pair correlation function 
of the SHSD fluid (solid symbols) , compared to MC (open sym- 
bols) and HNC calculations (solid lines) for the linear core sys- 
tem, at various densities and x = !• 
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Interestingly, in the limit x ~~ > 00 the SHSD algorithm 
reduces to hard-sphere (HS) molecular dynamics. In fact, 
if the density 4> is smaller than the freezing point for the 
HS system, the structure of the SHSD fluid approaches, as 
X increases, that of the HS fluid. For higher densities, if 
X is sufficiently high, crystallization is observed in SHSD, 
either to the usual hard-sphere crystals if <j> is lower than 
the close-packing density, or if not, to an unusual partially 
ordered state with multiple occupancy per site, typical of 
weakly repulsive potentials. 

An exact BBGKY-like hierarchy of Master equations for 
the s-particle distribution functions of the SHSD fluid is 
given in Ref. [17]. For the first equation of this BBGKY 
hierarchy, valid at low densities, we can neglect correlations 
other than pair ones and approximate /zfri, vi, r 2 , V2) = 
<?2( r i2)/i(ri, Vi)/(r2, V2). With this assumption we obtain 
a stochastic Enskog equation similar to a revised Enskog 
equation for hard spheres but with a smeared distribution 
of hard-sphere diameters, as studied in Ref. [18|. The 
Chapman-Enskog expansion carried out in Ref. [18j pro- 
duces the equation of state (EOS) p — PV/NkT, and ap- 
proximations to the self-diffusion coefficient £, the shear 77 
and bulk 7773 viscosities, and thermal conductivity n of the 
SHSD fluid. The expressions ultimately give the trans- 
port coefficients in terms of various integer moments of 
32(2;), x k = J Q x k g 2 (x)dx, specifically, p- 1 = 120x^3, 
C/Co = 0^/(480X^2), Vb/vo = 48<t> 2 xxi/ir 3/2 , and 



j]/rjo or k/k = 



ci 



(1 + c 2 4>XX3) 2 + c 3 r] B /vo, 



where £0 = Dy^kT/m, t?q = D 2 \J mkT and 
kD~ 2 ^/kT/m are natural units, and c\ — 5/48, c 2 — 24/5 
and C3 = 3/5 for ij, while ci = 25/64, C2 = 24/5 and 
C3 = 3/5 for k. 

The above formula for the pressure is exact and is equiv- 
alent to the virial theorem for the linear core potential, 
and thus thermodynamic consistency between g 2 (x) and 
p(<p) is guaranteed. In the inset in the top part of Fig. 
[2j we directly demonstrate the thermodynamic consistency 
of SHSD by comparing the compressibility calculated from 
the EOS, S c = (p + (fidp/dtfi)' 1 , to the structure factor 
at the origin Sq — S(lu — 0, k — 0). Furthermore, good 
agreement is found between the adiabatic speed of sound 
c 2 s — Sq 1 + 2p 2 /3 and the location of the Brilloin lines 
in the dynamic structure factor S(u>; k) for small k values. 
In Fig. [3J we also compare the theoretical predictions for 
77 utilizing the HNC approximation for g 2 to the ones di- 
rectly calculated from SHSD. Surprisingly, good agreement 
is found for the shear viscosity at all densities. The corre- 
sponding results for £ show significant (~ 25%) deviations 
for the self-diffusion coefficient at higher densities because 
of corrections due to higher-order correlations. 

As an illustration of the correct hydrodynamic behav- 
ior of the SHSD fluid and the significance of compressibil- 
ity, we study the velocity autocorrelation function (VACF) 
C(t) — (v x (0)v x (t)) for a single neutrally-buoyant hard 
sphere of mass m and radius R suspended in an SHSD fluid 




Figure 2: (Color online) Comparison between numerical results 
for SHSD at several collision frequencies (different symbols) with 
predictions based on the stochastic Enskog equation using the 
HNC g2(x) (solid lines). The low-density approximations are 
also indicated (dashed lines). (Top) Normalized equation of 
state. The inset compares the compressibility (pressure deriva- 
tive, dashed lines) to the structure factor at the origin S(k —* 0) 
(symbols) , measured using a direct Fourier transform of the par- 
ticle positions for small k and extrapolating to k = 0. (Bottom) 
The shear viscosity at high and low densities (inset), as mea- 
sured using an externally-forced Poiseuille flow. There are sig- 
nificant corrections (Knudsen regime) for large mean free paths 
(i.e., at low densities and low collision rates). 

of mass density p. This problem is relevant to the model- 
ing of polymer chains or (nano)colloids in solution, and led 
to the discovery of a long power-law tail in C (t) [2(| ■ 
Here the solvent-solvent particles interact as in SHSD. The 
solvent-solute interaction is treated as if the SHSD parti- 
cles are hard spheres of diameter D s , chosen to be some- 
what smaller than their interaction diameter with other 
solvent particles (specifically, we use D s = D/4) for com- 
putational efficiency reasons, using an event-driven algo- 
rithm [3j . Upon collision the relative velocity of the solvent 
particle is reversed in order to provide a no-slip condition 
at the surface of the suspended sphere 

HE! (slip bound- 
aries give qualitatively identical results). For comparison, 
an ideal solvent of comparable viscosity is also simulated. 

Theoretically, C(t) has been calculated from the lin- 
earized (compressible) fluctuating Navier-Stokes (NS) 
equations [l9j |. The results are analytically complex even 
in the Laplace domain, however, at short times an invis- 
cid compressible approximation applies. At large times the 
compressibility does not play a role and the incompress- 
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ible NS equations can be used to predict the long-time tail. 
At short times, t < t c = 2R/c s , the major effect of com- 
pressibility is that sound waves generated by the motion of 
the suspended particle carry away a fraction of the momen- 
tum, so that the VACF quickly decays from its initial value 
C(0) = kT/m to C(t c ) w kT/M, where M = m+2irR 3 p/3. 
At long times, t > t V i SC — ApRj^/Sij, the VACF decays as in 
an incompressible fluid, with an asymptotic power-law tail 
(kT/m)(8V3n)~ 1 (t/t V i SC )~ 3 / 2 , in disagreement with pre- 
dictions based on the Langevin equation (Brownian dy- 
namics), C(t) — (kT/m) exp (—6irRHilt/m). We have es- 
timated the effective (hydrodynamic) colloid radius Rh 
from numerical measurements of the Stokes friction force 
F = — 6ttRh rjv. 




Figure 3: (Color online) The velocity autocorrelation function 
for a neutrally buoyant hard sphere suspended in a non-ideal 
SHSD (x — 1) solvent at two densities (symbols), as well as an 
ideal I-DSMC solvent (cj> = 0.5, \ = 0.62, symbols), at short 
and long times (inset). For the more compressible (less viscous) 
fluids the long time tails are statistically measurable only up to 
t/tviac ~ 5. The theoretical predictions based on the inviscid, 
for short times, or incompressible, for long times, Navier-Stokes 
equations are also shown (lines). The diameter of the nano- 
colloidal particle is only 2.5D, although we have performed sim- 
ulations using larger spheres as well with very similar results. 
Since periodic boundary conditions were used we only show the 
tail up to about the time at which sound waves generated by its 
periodic images reach the particle, ti, = L/c B . 

In Fig. G2 numerical results for the VACF for an I-DSMC 
solvent and an SHSD solvent at two different densities are 
compared to the theoretical predictions. It is seen, as pre- 
dicted, that the compressibility or the sound speed c s , de- 
termines the early decay of the VACF. The exponent of the 
power-law decay at large times is also in agreement with the 
hydrodynamic predictions. The coefficient of the VACF tail 
agrees reasonably well with the hydrodynamic prediction 
for the less dense solvents, however, there is a significant 
deviation of the coefficient for the densest solvent, perhaps 
due to ordering of the fluid around the suspended sphere, 
not accounted for in continuum theory. 

In closing, we should point out that for reasonable val- 



ues of the collision frequency (\ ~ 1) and density (<f> ~ 1) 
the SHSD fluid is still relatively compressible compared to a 
dense liquid, c 2 s < 10. Indicative of this is that the diffusion 
coefficient is large relative to the viscosity as in typical DPD 
simulations, so that the Schmidt number S c — ?/(pC) 1 
is less than 10 instead of being on the order of 100-1000. 
Achieving higher c s or S c requires high collision rates (for 
example, \ ~ 10 4 is used in Ref. [13| | ) and appropriately 
smaller time steps to ensure that there is at most one col- 
lision per particle per time step, and thus a similar com- 
putational effort as in molecular dynamics. The advantage 
of SHSD is its simplicity, easy parallclization, and simpler 
coupling to continuum methods such as fluctuating hydro- 
dynamics (ioj . 
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